This notebook utilizes the data from Martin et al 2019 and processed using Seurat pipeline. Here we explore the data and generate figures to understand the CD ileal data quality and the cell type differences with inflammation. The data used here is availabe to download as rds file through our IBDGC DataCommons
Uncomment the installation line if the environment doesnot have those packages installed.
#install.packages(c("ggpubr", "Seurat", "ggplot2", "dplyr","tidyr","viridis"))
#devtools::install_github("rpolicastro/scProportionTest")
#Load the libraries
library(ggpubr)
library(Seurat)
library(ggplot2)
library(dplyr)
library(tidyr)
library(viridis)
setwd('/Users/mamtagiri/Downloads/CD ileal seurat')
test<-readRDS("/Users/mamtagiri/Downloads/CD ileal seurat/Immune.rds")
test$annotation<-test@active.ident
##### UMI plot for the clusters ####
meta<-test@meta.data
p1=ggplot(meta, aes(x=nCount_RNA)) + geom_histogram(binwidth=1000, fill="#69b3a2", color="#e9ecef", alpha=0.9)
p2=ggplot(meta, aes(x=annotation, y=nCount_RNA,fill=seurat_clusters))+geom_boxplot(show.legend = FALSE)+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1,size=15))
ggarrange(p1,p2, labels = c("UMI histogram","UMAP distribution"),common.legend = FALSE,ncol = 1, nrow = 2)
We explore how the cell contribution for each cell type is across the two conditions (inflamed, uninflamed)
##### Table with cell number per condition per cluster
table(test@active.ident,test@meta.data$group)
##
## inflamed uninflamed
## Epithelial 102 366
## Central Trm 5621 4244
## CD8 cytotoxic T cells 2665 4441
## IgA plasma cells 2955 3315
## gd T cells 2234 1025
## Naive B cells 2134 1070
## CD8+ Trm 951 1532
## Trm 1154 1318
## Mem. B cells 2716 1276
## Type 3 cytokine Trm 598 1223
## IgM plasma cells 674 932
## ILC1 830 651
## Tregs 1065 267
## DC2 614 687
## IgG plasma cells 2126 381
## Residential macrophages 535 582
## ILC3 472 634
## Fibroblasts 417 731
## moDCs 319 445
## ACKR1+ endothelial 486 219
## CD36+ endothelial 125 491
## Inflammatory macrophages 570 23
## DC1 335 153
## Activated fibroblasts 370 3
## Activated DCs 315 20
## Smooth muscle 132 128
## Glial cells 29 220
## Lymphatics 107 105
## Plasmablasts 133 59
## Pericytes 27 96
## UMAP of the 11 paired inflamed, uninflamed ileal single cell data ###
p1=DimPlot(test,reduction = "umap",label = TRUE,label.size = 6.5,repel = TRUE)+NoLegend()
p2=DimPlot(test,reduction = "umap",label = TRUE,label.size = 5,repel = TRUE,split.by = "group")+NoLegend()
ggarrange(p1,p2, labels = c("UMAP","UMAP by condition"),common.legend = FALSE,ncol = 1, nrow = 2)
#### Make a cell proportion bar graph ###
meta<-test@meta.data
counts <- as.data.frame(table(test@active.ident,test@meta.data$group))
ggplot(counts, aes(fill=Var2, y=Freq, x=Var1)) +
geom_bar(position="fill", stat="identity",alpha=0.8)+ scale_fill_brewer(palette="Set2") + xlab(" Clusters")+ ylab("Frequency") + labs(fill = "Condition")+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1,size=15))
library("scProportionTest")
prop_test <- sc_utils(test)
prop_test <- permutation_test(
prop_test, cluster_identity = "annotation",
sample_1 = "uninflamed", sample_2 = "inflamed",
sample_identity = "group"
)
permutation_plot(prop_test)
Visualize how consistent these inflammation related cell type changes are in each patient
counts <- as.data.frame(table(test@active.ident,test@meta.data$group,test$SampleID))
#Renaming cells under broader category
Plasma<-c("IgA plasma cells","IgG plasma cells","IgM plasma cells","Plasmablasts")
MNP<-c("Residential macrophages","Inflammatory macrophages","moDCs","pDCs","DC1","DC2","Activated DCs")
Tc<-c("Central Trm","CD8 cytotoxic T cells","gd T cells","CD8+ Trm","Trm","Type 3 cytokine Trm","Tregs")
Stromal<-c("Smooth muscle","Fibroblasts","Activated fibroblasts","Lymphatics","Pericytes","CD36+ endothelial","ACKR1+ endothelial","Glial cells")
B<-c("Mem. B cells","Naive B cells")
Mast<-c("Mast cells")
ILC<-c("ILC3","ILC1")
Epithelial<-c("Epithelial")
counts$id<-paste0(counts$Var3,"-",counts$Var2)
counts_new<-counts%>%mutate(CellGroup = case_when(
Var1%in%Plasma ~ "Plasma",
Var1%in%MNP ~ "MNP",
Var1%in%Tc ~ "Tcells",
Var1%in%Stromal ~ "Stromal",
Var1%in%B ~ "B cells",
Var1%in%Mast ~ "Mast",
Var1%in%Epithelial ~ "Epithelial",
Var1%in%ILC ~ "ILCs"
))
ggplot(counts_new, aes(fill=CellGroup, y=Freq, x=id)) + geom_bar(position="fill", stat="identity",alpha=0.8)+ scale_fill_brewer(palette="Set2") + xlab(" Clusters")+ ylab("Frequency") + labs(fill = "Condition")+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1,size=7))
### Overlay Cell type markers from Martin et al. 2019
featuresimmune<- c("CD3D", "CD2", "CD7", "TNFRSF17", "XBP1", "BANK1", "CD79B", "CD22", "MS4A1", "HLA-DRB1", "HLA-DQA1", "LYZ", "IL3RA", "IRF7", "GZMB", "LILRA4", "CLEC4C", "TPSAB1", "CMA1", "KIT", "FCER1A", "FCER1G", "KLRB1", "TYROBP", "PLVAP", "VWF", "COL3A1", "COL1A1", "ACTA2")
test<-ScaleData(test,features = featuresimmune)
DoHeatmap(subset(test,downsample=400),size=7,angle=90,lines.width=13,features = featuresimmune)+NoLegend() +scale_fill_viridis()+theme(axis.text.y = element_text(size = 15,colour = "black"))